package ij.hkust.QuickPALM;

import ij.*;
import ij.IJ.*;
import ij.measure.*;
import ij.plugin.*;
import ij.plugin.filter.*;
import ij.plugin.frame.*;
import ij.process.*;
import ij.gui.*;
import ij.measure.CurveFitter.*;
import java.awt.*;
import java.lang.Math;

import java.util.Random;

/**
 * 
 * @author linshan
 *
 */
public class Simulator implements PlugIn
{
	ImagePlus imp1;
	ImagePlus imp2;
	ImagePlus imp3;
	ImageProcessor ip1;
	ImageProcessor ip2;
	ImageProcessor ip3;
	ImageStack res;

	Random rd = new Random();

	MyDialogs dg = new MyDialogs();

	double [][] board;

	// double [] PoissonTable;

	public void run(String args)
	{
		IJ.register(Simulator.class);

		if (!dg.drawer())
			return;

		/*if (dg.sim_pssn)
			{
			loadPoissonTable(dg.sim_center);
			if (PoissonTable == null)
			{
			IJ.error("Can not load noise table");
			return;
			}
			}*/

		int owidth;
		int oheight;
		int width;
		int height;

		if (dg.sim_mode == dg.sim_modes[0])
		{
			owidth = dg.sim_clwidth;
			oheight = dg.sim_clheight;
			width = owidth * (int)(dg.sim_campsize*1000/(dg.sim_minpsize*dg.sim_magn));
			height = oheight * (int)(dg.sim_campsize*1000/(dg.sim_minpsize*dg.sim_magn));

			board = new double[owidth][oheight];

			res = new ImageStack(owidth, oheight);

			for (int index = 0; index < dg.drawer_nframe; index++)
			{
				ip3 = new ShortProcessor(owidth, oheight);

				for (int i = 0; i < owidth; i++)
				{
					for (int j = 0; j < oheight; j++)
					{
						board[i][j] = 0;
					}
				}

				SetNoise(owidth, oheight);

				if (dg.sim_obj == dg.sim_objs[0])
				{
					if (rd.nextDouble() < dg.sim_prob)
						SetSpot((int)(width/2), (int)(height/2), width, height, owidth, oheight, 1);
				}

				else if (dg.sim_obj == dg.sim_objs[1])
				{
					if (rd.nextDouble() < dg.sim_prob)
						SetSpot((int)(width/2), (int)(height/2), width, height, owidth, oheight, 1);
					if (rd.nextDouble() < dg.sim_prob)
						SetSpot((int)(width/2) + (int)(dg.sim_dis/dg.sim_minpsize), (int)(height/2), width, height, owidth, oheight, 1);
				}

				else if (dg.sim_obj == dg.sim_objs[2])
				{
					for (int m = (int)(width/4); m <= (int)(width*3/4); m++)
					{
						if (rd.nextDouble() < 0.5 * dg.drawer_intensity)
							SetSpot(m, (int)(height/2), width, height, owidth, oheight, 1);
						if (rd.nextDouble() < 0.5 * dg.drawer_intensity)
							SetSpot(m, (int)(height/2) + (int)(dg.sim_dis/dg.sim_minpsize), width, height, owidth, oheight, 1);
					}
				}

				else if (dg.sim_obj == dg.sim_objs[3])
				{
					int d = (int)(dg.sim_mindis/dg.sim_minpsize);
					for (int m = 0; (int)(width/5) + (int)(m*(m-1)*d/2) < (int)(width*4/5); m++)
					{
						for (int n = 0; (int)(height/5) + (int)(n*(n-1)*d/2) < (int)(height*4/5); n++)
						{
							if (rd.nextDouble() < dg.sim_prob)
								SetSpot((int)(width/5)+(int)(m*(m-1)*d/2), (int)(height/5)+(int)(n*(n-1)*d/2), width, height, owidth, oheight, 1);
						}
					}
				}

				else if (dg.sim_obj == dg.sim_objs[4])
				{
					int d = (int)(dg.sim_mindis/dg.sim_minpsize);
					for (int m = (int)(width/4); m <= (int)(width*3/4); m++)
					{
						for (int n = 0; (int)(height/5) + (int)(n*(n-1)*d/2) < (int)(height*4/5); n++)
						{
							if (rd.nextDouble() < 0.5 * dg.drawer_intensity)
								SetSpot(m, (int)(height/5)+(int)(n*(n-1)*d/2), width, height, owidth, oheight, 1);
						}
					}
				}

				Draw(ip3, owidth, oheight);
				res.addSlice(index + "", ip3, 0);
			}

			imp3 = new ImagePlus("Sample file", res);
			imp3.show();
		}

		else if (dg.sim_mode == dg.sim_modes[1])
		{
			if (dg.sim_open)
			{
				imp1 = IJ.openImage();
			}
			else
			{
				imp1 = IJ.getImage();
			}

			if (imp1 == null)
			{
				IJ.noImage();
				return;
			}

			ip1 = imp1.getProcessor().duplicate();

			width = imp1.getWidth();
			height = imp1.getHeight();

			if (dg.drawer_cusres)
			{
				owidth = (int)(width*dg.drawer_magn*dg.drawer_minpsize/(1000*dg.drawer_psize));
				oheight = (int)(height*dg.drawer_magn*dg.drawer_minpsize/(1000*dg.drawer_psize));
			}
			else
			{
				owidth = dg.drawer_owidth;
				oheight = dg.drawer_oheight;
			}

			board = new double[owidth][oheight];
			for (int i = 0; i < owidth; i++)
			{
				for (int j = 0; j < oheight; j++)
				{
					board[i][j] = 0;
				}
			}

			double magn = width / owidth;

			double prob[][] = new double[width][height];
			double max = ip1.getMax();

			for (int i = 0; i < width; i++)
			{
				for (int j = 0; j < height; j++)
				{
					prob[i][j] = ip1.getPixel(i, j) / max;
				}
			}

			if (dg.sim_climg)
			{
				ip2 = new ShortProcessor(owidth, oheight);

				for (int i = 0; i < owidth; i++)
				{
					for (int j = 0; j < oheight; j++)
					{
						board[i][j] = 0;
					}
				}

				SetNoise(owidth, oheight);

				for (int m = 0; m < width; m++)
				{
					for (int n = 0; n < height; n++)
					{
						if (prob[m][n] > 0.75)
							SetSpot(m, n, width, height, owidth, oheight, prob[m][n] * 0.1);

					}
				}

				Draw(ip2, owidth, oheight);
				imp2 = new ImagePlus("Sample file", ip2);
				imp2.show();
			}

			else
			{

				res = new ImageStack(owidth, oheight);

				for (int index = 0; index < dg.drawer_nframe; index++)
				{
					ip2 = new ShortProcessor(owidth, oheight);

					for (int i = 0; i < owidth; i++)
					{
						for (int j = 0; j < oheight; j++)
						{
							board[i][j] = 0;
						}
					}

					SetNoise(owidth, oheight);

					for (int m = 0; m < width; m++)
					{
						for (int n = 0; n < height; n++)
						{
							if (rd.nextDouble() < prob[m][n] * dg.drawer_intensity)
								SetSpot(m, n, width, height, owidth, oheight, 1);
						}
					}

					Draw(ip2, owidth, oheight);
					res.addSlice(index + "", ip2, 0);
				}

				imp2 = new ImagePlus("Sample file", res);
				imp2.show();
			}
		}
	}

	void SetNoise(int owidth, int oheight)
	{
		for (int i = 0; i < owidth; i++)
		{
			for (int j = 0; j < oheight; j++)
			{
				board[i][j] += (2 * Math.sqrt(dg.drawer_max) * dg.drawer_noise * rd.nextDouble());
			}
		}
	}

	void SetSpot(int x, int y, int width, int height, int owidth, int oheight, double fact)
	{
		int peak;

		if (dg.sim_climg)
		{
			peak = (int)(dg.drawer_max * fact);
		}
		else
			peak = (int)(dg.drawer_max - dg.drawer_peakdev * rd.nextDouble());

		if (peak < 0)
			peak = 0;

		int xstart;
		int ystart;
		int xend;
		int yend;

		double magn = width / owidth;

		double radius = magn * dg.drawer_radius;

		/*double [][] pvalues = new double[owidth][oheight];

			for (int i = 0; i < owidth; i++)
			{
			for (int j = 0; j < oheight; j++)
			{
			pvalues[i][j] = ip.get(i, j);
			}
			}*/

		if (dg.drawer_mode == dg.drawer_modes[0])
		{
			xstart = (x - (int)(1.5*radius) > 0)?(x - (int)(1.5*radius)):0;
			ystart = (y - (int)(1.5*radius) > 0)?(y - (int)(1.5*radius)):0;
			xend = (x + (int)(1.5*radius) < width)?(x + (int)(1.5*radius)):width - 1;
			yend = (y + (int)(1.5*radius) < height)?(y + (int)(1.5*radius)):height - 1;

			for (int i = xstart; i <= xend; i++)
			{
				for (int j = ystart; j <= yend; j++)
				{
					double r = Math.hypot(i - x, j - y);
					double add = peak*Math.exp(-2*r*r/(radius*radius));
					board[(int)(i/magn)][(int)(j/magn)] += add/(magn*magn);
				}
			}

			/*for (int i = 0; i < owidth; i++)
				{
				for (int j = 0; j < oheight; j++)
				{
				int old = ip.get(i, j);
				int add = (int)(Math.round(pvalues[i][j]/(magn*magn)));
				ip.set(i, j, old + add);
				}
				}*/

		}
	}

	void Draw(ImageProcessor ip, int owidth, int oheight)
	{
		if (ip == null)
			ip = new ShortProcessor(owidth, oheight);

		for (int i = 0; i < owidth; i++)
		{
			for (int j = 0; j < oheight; j++)
			{
				ip.set(i, j, (int)board[i][j]);
			}
		}
	}

	/*double pnoise()
		{
		double r = rd.nextDouble();
		if (PoissonTable != null)
		{
		for (int i = 0; i < 100; i++)
		{
		if (r > PoissonTable[i]/PoissonTable[99] && r < PoissonTable[i+1]/PoissonTable[99])
		return i/100;
		}
		}
		return r;
		}

		void loadPoissonTable(double center)
		{
		if (center <= 0 || center >= 1)
		{
		IJ.error("Wrong parameter input, return default value(0.5)");
		PoissonTable = null;
		return;
		}

		int mean = (int)(center * 100);

		PoissonTable = new double[100];
		double sum;

		for (int i = 0; i < 100; i++)
		{
		sum = 0;
		for (int j = 0; j < i; j++)
		sum += poisson(j+1, mean);
		PoissonTable[i] = sum;
		}
		}  

		double poisson(int i, int mean)
		{
		return exponent(mean, i)/(Math.exp(mean)*factorial(i));
		}

		int factorial(int i)
		{
		int product = 1;
		for (int j = 1; j <= i; j++)
		product *= j;
		return product;
		}

		int exponent(int i, int j)
		{
		int product = 1;
		for (int k = 0; k < j; k++)
		product *= i;
		return product;
		}*/
}


